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Abstract 

Finite difference methods for solving problems of time-harmonic 
acoustics are developed and analyzed. Multi-dimensional inhomoge- 
neous problems with variable, possibly discontinuous, coefficients are 
considered, accounting for the effects of employing non-uniform grids. 
A weighted-average representation is less sensitive to transition in 
wave resolution (due to variable wave numbers or non-uniform grids) 
than the standard pointwise representation. Further enhancement in 
method performance is obtained by basing the stencils on generaliza- 
tions of Pade approximation, or generalized definitions of the deriva- 
tive, reducing spurious dispersion, anisotropy and reflection, and by 
improving the representation of source terms. The resulting schemes 
have fourth-order accurate local truncation error on uniform grids and 
third order in the non-uniform case. Guidelines for discretization per- 
taining to grid orientation and resolution are presented. 
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1 Introduction 


Boundary- value problems governed by the Helmholtz equation model propa- 
gation and evanescence of time-harmonic waves, describing a variety of physi- 
cal phenomena, including acoustic, elastic and electromagnetic waves. When 
the wavelength is of the same order as characteristic length scales asymptotic 
methods usually cannot be employed and standard numerical methods such 
as boundary element, finite element or finite difference methods are often 
required. 

Finite difference methods are a prevalent computational technique that 
applies to variable coefficient as well as nonlinear problems. A general frame- 
work for deriving higher-order finite difference schemes was proposed by 
Lynch and Rice for ordinary differential equations [1] and elliptic partial 
differential equations [2], and applied to the Helmholtz equation [3]. The 
coefficients of the stencil in this method are computed by solving a local sys- 
tem of equations so that it is exact on a given space of polynomials. This is 
repeated at every gird point at which the solution is unknown, contributing 
to the cost of computation. Accuracy is further enhanced by judiciously se- 
lecting the points at which source terms are evaluated and computing their 
coefficients in the same way. 

A family of finite difference schemes with improved representation of a 
range of wave numbers is presented and analyzed in [4], Tam and Webb [5] 
optimize the dispersion properties of higher-order finite difference discretiza- 
tion of the linearized Euler equations. In this approach the order of accuracy 
of a stencil is allowed to drop, freeing a parameter that is then designed to 
improve dispersion response. 

Finite difference equations can be obtained by replacing the limit that 
defines the derivative with a finite counterpart. The familiar definition of 
the derivative may be generalized by introducing a resolution-dependent pa- 
rameter leading to improved performance of the discrete methods. As long as 
the limit behavior is unaltered consistency of the approximation is retained. 
This idea was introduced by Mickens and employed as a means of generat- 
ing stable finite difference schemes on uniform grids for several differential 
equations in one spatial dimension ([6, 7] and references therein). Similar 
discrete equations are obtained by new classes of finite element methods for 
a variety of applications, including wave propagation (e.g., [8] and references 
therein). It should be noted that analysis of the computation of waves [9] 
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indicates that accuracy depends not only on the product of wave number 
and grid size (related to resolution), but also on product of other powers of 
these quantities. 

In this work we apply generalizations of the definition of the derivative and 
of Pade approximation, to finite difference stencils for the Helmholtz equa- 
tion in order to obtain improved dispersion behavior. Contrary to HODIE 
methods, the coefficients are obtained explicitly. Multi-dimensional inhomo- 
geneous problems with variable, possibly discontinuous, coefficients are con- 
sidered, accounting for the effects of employing non-uniform grids. Several 
finite difference stencils in one and two dimensions are presented in Sec. 2. 
The analysis of the numerical methods gradually builds up to the general 
case. Performance of the various formulations for homogeneous problems 
with constant coefficients on uniform grids is examined by dispersion anal- 
ysis in Sec. 3. This tool is employed to define the resolution-dependent 
parameter for improved performance. In Sec. 4 the effect of the direction 
of wave propagation relative to grid lines is accounted for. The effects of 
non-uniform grids and discontinuities in physical properties are investigated 
in Sec. 5. Standard finite difference methods are often not well-suited for 
interface problems (see, e.g., [10, pp. 17-21]). However, appropriate repre- 
sentation preserves the order of accuracy of the continuous-coefficient, and 
even constant-coefficient case. (Issues related to curved interfaces, as well as 
curved domain boundaries are not treated herein.) The results of these anal- 
yses are corroborated by means of local truncation error analysis in Sec. 6, 
accounting also for the effects of source terms. Numerical testing of these 
stencils will be performed in future work. 


2 Discrete Formulations 

The Helmholtz equation is 


&<f> + k 2 (j> + / = 0 ( 1 ) 

where k = lj/cq is the wave number, oj is the angular frequency and Cq is the 
speed of sound, and / is a given source distribution. Although not explicitly 
addressed in the following, the case of k 2 < 0 which corresponds to evanescent 
waves or singular diffusion problems is also covered. An inhomogeneous 
medium is represented by spatial variations in k[x ). 
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2.1 One dimension 


Consider a uniform grid of size h with points at x 3 — jh. A typical start- 
ing point is based on the standard finite difference stencil for the second 
derivative 


& XX 


<f>j + 1 “ 2 <t>j + <t>j - 1 
h 2 


( 2 ) 


and pointwise (PT) representation of undifferentiated terms 


D xx (f)j + k 2 (f>j + /j — 0 


( 3 ) 


where <j> 3 is the discrete solution at point j and fj = /(xj). On a non-uniform 
grid this generalizes to 


/ (j>j 

V h + 


4>j - 4>J-1 \ / + h 

h^~ 


+ + fj — 0 


( 4 ) 


where h and h+ are the grid size before and after point j , respectively. For 
a discontinuity in physical properties at point j the stencil becomes 


( $3 

V 



'*+ + *- , {k+) 2 h+ + (k-) 2 h- , , , 

2 + A+ + /i- 


0 (5) 


where k~ and are the wave numbers before and after point j, respec- 
tively. These may be also considered as piecewise-constant approximations 
of variable coefficients. 

The undifferentiated terms may be represented by a weighted average 
(wa) suggested by linear finite elements (with piecewise linear approximation 
of the source distribution, see, e.g., [11, pp. 45-46]) 

+ e 6+.+4fr + + /,■+■ + Vi + Si-) = o (6) 

6 o 

This scheme has the same asymptotic behavior as the pointwise representa- 
tion, but improvement in performance on coarse grids is evident (see Sec. 3). 
For variable coefficients this becomes 

r, 1 , [k 2 (j))j+\ + \{k 2 (i>)j + (k 2 <f>)j-i fj+i + Afj + fj-1 m 

L>xx<t>] + g 1- q u \ < ) 

where {k 2 (f))j = k 2 (xj)4>j. 
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On a non-uniform grid the appropriate weighting is 


P 

y 


y+1 - <pj / h + + h: 


h+ 


h- 


2 


+ 


(</>j+i + 2 (j>j) + — (2 <j)j -I- + 


K h+ + h~ 
1 ( h + 


h+ + h- 


3 \h+ + h 


r(/i+i + 2/j) + 


ft+ + ft 


r (2/j + /j-i) 


0 (8) 


Superior performance on non-uniform grids (see Sec. 5) is attained with no 
increase in the number of points in the stencil. For a discontinuity in physical 
properties at point j the stencil becomes 


1 

3 


/ <j ) J+ , - (f>j - ^ / ft+ + ft’ 

\ h+ h- 


+ 


[k * Yh * z (^ + 2*) + £££(2*; + + 


^ ft+ + ft 
1 / ft + 


3 \ft+ + ft- 


-(/j+i + 2/,) + 


ft+ + h- 
h- 


ft+ + ft 


r(2 fj + /j-i) 


0 


( 9 ) 


Again, this may also be considered as a piecewise-constant approximation of 
the case of variable coefficients. 

Performance of finite difference schemes for the Helmholtz equation may 
be enhanced by basing the stencils on more general definitions of the deriva- 
tive 

d<f> <j>(x + ft) — <f>(x ) 

— = lim — —t— —Z- 

dx h->o fi(kh) h 

where, for consistency 

lim $ = 1 (11) 

kh-* 0 v ’ 



This definition depends on kh, an indication of wave resolution by the grid. 
For the Laplacian this reduces to the standard definition for grids of any size. 

This generalization of the derivative definition may be applied for either 
the first or second derivatives, or to both. On uniform grids all are equivalent. 
Since the parameter depends on the grid size it is applicable to non-uniform 
grids as well. Superior performance on non-uniform grids is obtained by 
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applying this concept to the second derivative alone (see Sec. 5). For the 
uniform case this reduces to 

- 2(f>j + d'j-i , 2 <h± ] + 4< ftj + < tlzl | h±l + 4 -^ + -fc- 1 - Q (12) 
fih 2 6 6 

The resolution-dependent parameter fi is defined to improve method per- 
formance. For example, the parameter may be defined to eliminate numerical 
dispersion 

a 6 1 ~ COs(fcft) /jgx 

( kh ) 2 2 + cos (kh) 

so that in simplified settings the phase is exact (EP), resulting in no trunca- 
tion error under some circumstances. This definition satisfies the consistency 
requirement (11). In such cases the representation of source terms, which is 
exact for piecewise linear variation, is no longer sufficiently accurate. A 
modification of the representation of source terms that does not degrade the 
higher-order accuracy of such schemes, similar to that employed by HODIF 
methods [2], is 

4 > j+ 1 — 2,<f>j + 4>j - 1 ,2 &}+] T + 4 >)- 1 /j+i/2 T /; T fj- 1/2 _ (14) 

Jh 2 6 3 

(suggested by linear finite elements with piecewise quadratic approximation 
of the source distribution) where /j+1/2 is the load term evaluated at the 
midpoint. For a piecewise linear source distribution this is identical to (12). 

One possibility of the parameter 


0 = 


12 

12 - (kh) 2 


(15) 


yields high-order representation 

Hji $ j - f" k 


2 <f>j+ 1 + 10 4>j + <j>j- 1 /y+1/2 + fa + /j-1/2 _ 

+ 3 


12 


= 0 


(16) 


(ho). This stencil (without the modification in the representation of source 
terms) may also be derived by employing Pade approximation 


k) XX 

1 T h 2 /\2 D x 


-<f>j + k 2 <j)j + /j — 0 


(17) 



(see, e.g., [12, p. 538]). 

This concept, in its original form, which may be viewed as an average of 
the pointwise and weighted-average representations of the undifferentiated 
term [13], provides high-order performance on uniform grids, but severely 
degrades in the non-uniform case. However, an appropriate generalization to 
non-uniform grids, based on the concept of generalizing the derivative defi- 
nition, leads to improved performance in the general case as well (see Secs. 5 
and 6). Allowing discontinuities in physical coefficients and accounting for 
non-uniform grids the proposed scheme is 


1 

3 

1 

3 


/ ^j + 1 - 4*j _ <t> 3 - <f>j- 1 
\ h+ h- 

( (*+) 2 0+A+ ,, , x . (k-y/3-h- 


' 0+ h++p~h- 


+ 


-{2<f)j + I + 


{/3+h+ + 0-h- ( ^ +1 + + 0+h+ + fi-h- 

\J+h+ + p-h-^ i+1/7 + ^ + ]}+h+ + + 2 ^ -1/2 ) 


where /?* = /?(A; :t /t ± ). 


(18) 

0 


2.2 Two dimensions 

Consider a two-dimensional uniform grid of size h with points at Xi = ih 
and y } = jh. For simplicity we consider the homogeneous case. A typical 
starting point is the five-point representation 


~ 2<f>jj ~h <f>i - i 4*i,j+i 


D XX 4*1,3 T Dyy4>i j T k 4*1,3 

24*1,3 + 4*i,j - 1 


h 2 


h 2 


+ k*4> lt] = 0 (19) 


which is the two-dimensional analog of (3). Non-uniform grids and material 
discontinuities may be accounted for by generalizations of (4) and (5). 

The two-dimensional counterpart of the idea that leads to (6), obtained 
by employing bilinear finite elements, is a nine-point representation 

^i+i j-i + 24 >j+\,j + d’t+ij+i — 8<ft,j -f fa-ij-i + 2<j*i- i,j + 4*j- i,j-n 

h 2 

^i-i, j+i + 24. i.i+i T d>;+i,j+i — 8 4>i,j T 4*i-\,j-\ + 2 4*i,j-\ + 4i+\,3-\ 

h 2 ’ + 
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— + 4 4- + <&+ij-i + 160ij + (20) 

o 

+ 4<^tj_i + 4^,-ij + ^t-ij-i) = 0 

leading to a significant reduction in spurious anisotropy (see Sec. 4). HODIE 
methods [2] also employ nine-point stencils in two dimensions. The band- 
width of the ensuing linear algebra problem is typically slightly larger but 
the difference in the cost of computation is insignificant. 

Performance can again be improved by substituting flh for h as in the one- 
dimensional case (14), based on the same definitions (13) and (15), although 
the methods are higher order only for propagation along grid lines. 

In order to maintain higher-order performance on uniform grids in two 
dimensions in all directions of propagation, the Pade approximation concept 
is employed. The two-dimensional counterpart of (17) is 


1 + 


D x 

h? 


12 Dxx 


4>i J + 


D 


yy 


1 + 


— D 

12 w 


+ k 2 — 0 


This may be generalized to 

^1 T Y 2 D xx <f>ij + ^1 T Y 2 Dyy ( Pi,] f" 

k 2 ^ 1 + —(D XX + Dyy ) + 7 D XX Dyy^j <f>ij = 0 


( 21 ) 


( 22 ) 


where 7 is selected to further improve properties in directions other than 
along grid lines, without effecting dispersion along grid lines and without 
degrading higher-order behavior in all directions. 

The standard Pade approximation is obtained by selecting 7 = 1 which 
yields the scheme 

+ Dyy H — D XX Dyy ^ (f > , , j T 

+ 100<fcj+ (23) 

<t>i~ i,j+i + 10^>, ,j_i + 10^ t _ij + <£,-ij-i) = 0 
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where 


^Dxx + Dyy + — D XX Dyy^j 4, J 

4+i, j-i + 84+1,7 + 4+ij+ i - 204 , j + 4-ij-i + 84-i , j + 4-ij+i 

h 2 

4-1, j+i + 84,7+1 + 4+ij+i - 204, j + 4-i, j-i + 84, ,-1 + 4+i, i-i 


+ 


(24) 


Neglecting higher-order terms in the Pade approximation by selecting 7 = 0 
leads to the slightly simplified stencil presented in [12, p. 542] 

yDxX + Dyy + — D xx Dyy^j 4 ,7 + 

* 2 (4+i,7 + 4j+i + 84 (j + 4 >J _ 1 + 4— 1 ,j ) = 0 (25) 

The computational cost is essentially unaffected since the bandwidth of the 
algebraic equations is identical. Another alternative presented in [12, p. 542] 
is 


yD XX T Dyy + g D xx Dyy^j 4jT 

k 2 

■g" (4+ij+i + 44+17 + 44,7+1 + 4+i, j-i + 524,7+ (26) 

4-i, 7+1 + 44,7-1 + 44_,j + 4-u-i) = 0 

obtained by selecting 7 = 2. 

Other values for 7 lead to other alternatives. In Sec. 4 it is seen that 
selecting 7 = 2/5, which leads to the stencil 

i^Dx X + Dyy + — D XX Dyy ^ 4 ,j + 

k 2 

2 q ( 4+i, 7+1 + 284+ij ~h 284,7+1 + 4+i, j-i + 2444,7+ (27) 

4-i, 7+i + 284,7-1 + 284-1,7 + 4-1, i-i ) = 0 

minimizes dispersion along the diagonals. On the other hand, reducing sen- 
sitivity of the scheme to direction of propagation is attained by the choice of 
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7 = 14/5, which yields the stencil 


4 “ D yy ~\~ ^ D XX DyyJ 

k 2 

^(7<^i+i j+i + + 16</>^j+i + 7 <^ 2 - 4-1 i + 268</> 1J + (28) 

7<f>i- ij+i + j-i ) = 0 

All these alternatives reduce to (ho) in one dimension. Thus the dispersion 
analysis for (ho) in vSec. 3 describes the dispersion of all alternatives along 
grid lines. In Sec. 4 the performance of various alternatives in other directions 
is compared. 


3 Spurious Dispersion 


A homogeneous, isotropic continuum is non-dispersive. 
phase velocity 

to 


c p 



and energy propagates at the group velocity 


dto 

~dk: 


— co 


Waves travel at a 

(29) 

(30) 


and so both are identical. 

For the discrete representation this is usually no longer the case. Each 
stencil can be characterized by an approximate wave number k h = k h (kh ), 
which depends on wave resolution and thus accounts for numerical dispersion. 
The phase velocity in the discrete case is 


< 4 - 


u 

k h 


TT ~ 77 


k h 


and the numerical group velocity is 


c„ = 


du> 

W 

dio dk 


Ok dk h 
(dk^ _1 
\~dk 


Co 


(31) 


(32) 
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On a uniform grid in one dimension a numerical solution in the form of 
a plane wave is 

<j>j = exp (ik k h) J (33) 

PT For point j the pointwise representation (3) of the plane wave solution 
yields 

0 = exp(ik k h) — ^2 — (kh) 2 ^ + 1 / exp(ik h h ) 

= 2 cos(k h h) — ^2 — (kh) 2 ^ (34) 

leading to the dispersion relation 

k h h = arccos (\ — (kh) 2 / 2^ (35) 

In one dimension the number of grid points in a wave is 

G = 2 7r /(kh) (36) 

The discrete solution represents propagation in the range kh < 2 which 
corresponds to a limit of approximately three grid points per wave- 
length. Within this range the numerical phase velocity is thus 

Cp/co = kh/ arccos ^1 — (kh) 2 / 2^ (37) 

and the numerical group velocity is 

C J/co = y/l -(khy/ 4 (38) 

Both are slower than the speed of sound in the material cq. 


WA S imilarly, for the weighted-average representation (6) the dispersion 
relation is 


k h h = 


arccos 


/ l-(M) 2 /3 \ 
\i + {khy/6j 


(39) 


representing propagation in the range kh < y/\2, a limit of approxi- 
mately two grid points per wavelength. Within this range the numerical 
phase velocity is obtained directly from the dispersion relation and the 
numerical group velocity is 


cj/co = \f\ - (kh) 2 j\2 (l + (kh) 2 / 6) (40) 

Both are faster than the speed of sound in the material. 
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EP The resolution-dependent parameter (3 may be defined so that discrete 
representations are non-dispersive (13) as is the case for the continuum. 
In one dimension this formulation has zero local truncation error for 
the homogeneous, constant coefficient case on uniform grids, and the 
phase and group velocities are exact. Careful generalization leads to 
improved performance on general configurations. 


HO The higher-order representation (15) is an approximation of the exact 
phase definition (13) on uniform grids. The resulting higher-order dis- 
persion relation 


k h h = 


arccos 


/I — 5 (kh) 2 /12\ 

v 1 + (khy/\ 2 ) 


(41) 


is a (1,1) Pade approximation, representing propagation in the range 
kh < \/6> a limit of approximately 2 1/2 grid points per wavelength. 
Within this range the numerical phase velocity is again obtained di- 
rectly from the dispersion relation and the numerical group velocity 
is 

cj/co = \A - W )' 76 (l + (kh) 2 / 12) (42) 

Both are slower than the speed of sound in the material. The power 
series expansion of the dispersion relation 


k h h 


+ <«) t 


480 


12096 


> kh 


demonstrates the higher-order nature of this representation. 


Dispersion of the various formulations is plotted in Fig. 1. Note that the 
region of primary interest is G > 4, a resolution of at least four grid points 
per wavelength. Within this region the errors in the pointwise and weighted- 
average representations are similar, and the asymptotic behavior is the same 
(see Sec. 6). However, even approaching the limit of resolution, and certainly 
beyond it, the weighted-average performance is superior. For example, at 
the limit of resolution ( G = 4) there is a 38% error in the group velocity for 
the pointwise representation, whereas in the weighted-average representation 
the error is only 26%. The higher-order method offers significantly superior 
representation, an error of only 7% in group velocity at G = 4, and the 
exact-phase method provides dispersion-free solutions. 








4 Spurious Anisotropy 

On the uniform grid in two dimensions a numerical solution in the form of a 
plane wave oriented at angle 9 to the grid lines is 

<f>ij = exp (ik h h cos 0) 1 exp (ik h h sin 0) J (44) 

PT For point i , j the pointwise representation (19) of the plane wave solution 
yields 

0 = exp {ik h h sin 0 ) + exp (ik h h cos 9) — (4 — (A*/?) 2 ) + 

1 / exp(ik h h cos 9) + 1/ exp(ik h h sin 9) 

= 2 (cos(A ,/j 7r cos 0) + cos(A: /l /t sin 0)) — ^4 — (kh) 2 ^ (45) 

leading to various dispersion relations, depending on the angle of orien- 
tation 9. When the wave is aligned with the grid (e.g., 0 — 0) this leads 
to the one-dimensional dispersion relation (35). The other extreme case 
occurs when the wave is oriented in the direction of cell diagonals (e.g., 

9 = tt/4) 

k h h — \/2arccos (^1 — (kh) 2 /ij (46) 

The discrete solution represents propagation in the range kh < y/8. 
Within this range the numerical phase velocity is again obtained di- 
rectly from the dispersion relation and the numerical group velocity 

cj/c o = V 1 - (*M7 8 (47) 

Both are slower than the speed of sound in the material. It is interesting 
to note that the pointwise representation is more dispersive when waves 
are oriented along the grid. 

WA Similarly, for the weighted-average representation (21) dispersion rela- 
tions are obtained from 

4 (ft — ( kh ) 2 ) — (l2 + (fc/i) 2 ) cos(k k h cos 9) co$(k k h sin 9)— (48) 

2 (3 + (A:/*) 2 ) {cos{k h h cos 0) + cos(k h h sin 0)) = 0 
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When the wave is aligned with the grid this leads to the known one- 
dimensional dispersion relation (39) and for waves parallel to cell diag- 
onals 


k h h = \/ 2 arccos 


( 1 ~ (kh) 2 /6 \ 
\i + (khy/i2) 


m 


representing propagation in the range kh < \/24. Within this range 
the numerical phase velocity is obtained directly from the dispersion 
relation and the numerical group velocity is 


cj/co = 0 -(^) 2 /24 (l + (kh) 2 / 12 ) (50) 

Both are faster than the speed of sound in the material. The weighted- 
average representation is also more dispersive for waves that are ori- 
ented with the grid. 

HO Dispersion relations for the higher-order representation with / 3 defined 
in (15) at various angles of orientation may be found in similar fashion. 
For waves aligned with cell diagonals the relation is identical to that of 
the pointwise representation at this angle (46), and the same holds for 
the wave velocities. This representation is thus higher order only for 
waves oriented along grid lines. 

Representations that are truly higher order in all directions of propa- 
gation are based on ( 22 ). For waves aligned with the grid this leads to 
the higher-order dispersion relation (41) and along cell diagonals the 
relation is 


k k h — \/ 2 arccos 


'6^/1 + (1 — 7)(M) 4 / 144 - (4 + (6 - 7 )(M) 2 /12)' 


2 + ~f(kh) 2 /\2 


By examining the power series expansion of this relation 


(51) 


k h h « kh + (5 7 - 2 )-^—— f 

v ' ’ 5760 96768 


(52) 


it is clear that the value of 7 = 2/5 minimizes dispersion in the direction 
of cell diagonals. On the other hand, 7 = 14/5 minimizes the difference 
between the dispersion along grid lines and in the direction of cell 
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diagonals, as seen by comparison to the power series expansion of the 
higher-order dispersion relation along grid lines (43). Dispersion in the 
direction of cell diagonals for various values of 7 is plotted in Fig. 2. As 
expected, the stencil with 7 = 2/5 is essentially 11011 -dispersive in the 
region of primary interest with a resolution of at least four grid points 
per wavelength. 

The ratio between numerical dispersion along grid lines and along cell 
diagonals is shown in Fig. 3. By design, the stencil with 7 = 14/5 is the 
least anisotropic in the range of at least four grid points per wavelength. 
This is corroborated by the polar plots in Fig. 4, showing the variation 
in phase velocity with angle of orientation for various resolutions. Note 
that the figure shows (Cp/co ) 4 to accentuate deviations from the exact 
value of unity. As expected, all the schemes perform identically along 
grid lines, but behavior in other directions is determined by the choice of 
7 . Differences among the various cases become more pronounced with 
reduced resolution, but in general are not extreme. The two schemes 
that stand out are indeed 7 = 2/5 which minimizes dispersion along 
diagonals, and hence overall, and 7 = 14/5 which reduces anisotropy. 

EP The case that is non-dispersive in one dimension (13) may be treated 
similarly. In this case there is no dispersion for waves aligned with 
the grid and the dispersion relation for waves in the direction of cell 
diagonals is 

k h h = \/ 2 arccos 

The numerical phase velocity is again obtained directly from the dis- 
persion relation and the numerical group velocity is 


+ 2 cos( kh)' 
5 + cos (kh) 


(53) 


Cg 7 + 5 cos (kh) fb + cos(kh)\ 
Cq ^ 1 T cos (kh) y \/2 1 6 / 


(54) 


This representation is obviously less dispersive for waves that are ori- 
ented with the grid. 


Dispersion in the direction of cell diagonals of the various formulations 
is plotted in Fig. 5 . Recall that the region of primary interest is G > 4, a 
resolution of at least four grid points per wavelength. Dispersion properties of 







Phase velc 




Figure 3: The ratio of numerical dispersion along grid lines and along cell 
diagonals for higher-order discrete formulations based on generalized Pade 
approximation. 






each scheme at arbitrary orientations are bounded on one hand by dispersion 
along grid lines shown in Fig. 1 , and on the other hand by dispersion along 
cell diagonals shown in Fig. 5. Performance of the point-wise and weighted- 
average schemes improves the farther the orientation of propagation is from 
the direction of grid lines. The same holds for the higher-order schemes of 
interest with 7 = 2/5 and 7 = 14/5. For the (EP) method the opposite 
occurs, so that performance of this scheme is vastly superior along grid lines. 
This scheme is higher order only along grid lines, but it still maintains a high 
degree of phase accuracy in all orientations. 

The resolution-dependent parameter ft may be defined so that the nu- 
merical representation is non-dispersive for waves at any given angle of ori- 
entation. For example 

= _12_ 1 - COS (y/2kh/2) 

1 (khy 2 + cos (y/2kh/2) 

eliminates dispersion of waves along cell diagonals. Similar performance was 
attained in the context of finite element methods [14]. In general, how- 
ever, the direction of wave propagation is not known in advance and there 
is a concern that defining (5 for any orientation other than along grid lines 
may degrade performance on non-uniform grids, as discussed in the follow- 
ing section. Furthermore, grids should be aligned with dominant directions 
of propagation to the extent possible. For these reasons it is preferred to 
maintain dispersion-free discrete solutions along grid lines. 

Numerical dispersion is thus sensitive to the orientation of wave propaga- 
tion. The two extreme cases are along grid lines shown in Fig. 1 , and along 
cell diagonals shown in Fig. 5. The largest change in dispersion properties 
possible is thus the ratio between the two, shown in Fig. 6 . Recall that the 
region of primary interest is G > 4 , a resolution of at least four grid points per 
wavelength. For highly resolved phenomena the performance of all schemes 
is similar and quite good. As wave resolution is reduced only the higher- 
order schemes (with values of 7 shown) retain a low level of anisotropy. Of 
the other schemes, approaching the limit of resolution and certainly beyond 
it, the pointwise method is clearly more sensitive to direction of propaga- 
tion. For example, at the limit of resolution {G — 4), there is 8 % anisotropy 
in the pointwise representation of phase velocity, whereas the anisotropy of 
other methods is at most about half that value. This becomes even more 
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pronounced in group velocity. 

Figure 7 shows the variation in phase velocity with angle of orientation 
of different schemes at various resolutions. For presentation purposes the 
figure shows (c£/c 0 ) 4 . Note that the plot for (pt) does not include the case 
of G = 3 since this scheme no longer represents propagation at this low 
resolution, which, in any event, is outside the region of primary interest of 
G > 4 . It is clear from these plots that the numerical phase velocity is less 
than the speed of sound in the material in all cases shown except for (wa). 
Close examination of Fig. 4 indicates that this is true of higher-order methods 
only with 7 > 2/5. With the exception of (ep), all the schemes considered 
exhibit superior dispersion behavior along cell diagonals. This would not hold 
for higher-order methods with 7 > 14/5, but there is no apparent motivation 
to pursue such methods in the first place. As mentioned, employing (55) 
eliminates dispersion along cell diagonals, leading to a version of (ep) with 
superior dispersion behavior along diagonals that is similar to other schemes 
in this regard. 

Overall, high wave resolution or higher-order methods are required if 
anisotropy is a concern. Of the methods that are not high order, on grids with 
lower resolution, the weighted average representation and its enhancements 
are much less anisotropic than the standard point wise representation. 


5 Spurious Reflection and Transmission 

Reflected and transmitted waves are generated by incident waves at disconti- 
nuities in physical properties. Numerical dispersion of discrete formulations 
gives rise to incorrect representation of these phenomena at transitions in 
wave resolution. 

5.1 Grid transition 

In a homogeneous material no reflection should occur. However, changes in 
grid size alter wave resolution giving rise to spurious reflection and transmis- 
sion due to numerical dispersion, phenomena that may be characterized in 
a manner similar to that of waves at discontinuities in physical coefficients. 
These phenomena are well known [15], have been carefully analyzed [16] and 
numerically demonstrated [17]. 
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Consider a one-dimensional configuration discretized by a uniform grid 
of size h~ to the left of the origin and a uniform grid of size to its right. 
Due to the transition in grid size at the origin, an incident plane wave of unit 
magnitude traveling in the positive direction with discrete values given in 
(33) for j < 0 generates spurious transmission such that <^ 0 7 ^ 1 and spurious 
reflection of magnitude (j>o — 1 . The numerical solution is thus 




exp (ik n h ) J + (<f>o — l) / exp(ik h h ) J , j< 0 
<f> 0 exp(ik h * h + y , j > 0 


(56) 


where the dispersion error is represented by the numerical wave numbers 
k h = k h (kh ± ), and the transmission error is represented by <^ 0 - In particular 


4>\ = 4>o exp(ik h+ h + ) 

d»-i = <f>oexp(ik h h~) — 2is\n(k h h.~) 


(57) 

(58) 


PT The pointwise representation (4) of this solution at the origin yields 
2 


0 = 


h+ + Ir 

2 

h+ + h- 


4>\ ~ (1 - ( kh + y/2)<j>o - (1 - (klr ) 2 /2)<j>o' 

h+ + h~ 

(f>o exp(ik h+ h + ) — cos(k h+ h + )<f) 0 


h+ 


+ 


(f) 0 exp(ik h h ) — '2is\n(k h h ) — cos (k h h~)<j > 0 ' 


h + + fi- 


ll ~ 

"sin (k h+ h + ) sin (k h h~) 

h+ + h~ 


< f > 0 ~ 2- 


. sin(k /l h )' 


i(59) 


where the second line was obtained by the dispersion relation for the 
pointwise representation (35). Thus 


0o — 


2 yj\ - (kh~) 2 /4 


yj\ - ( kh^yji + y/l - (kh+y / 4 


(60) 


which is valid in the range of resolution in which the pointwise formu- 
lation represents propagation (along grid lines). 
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WA Similarly, the weighted-average representation (8) of the solution (56) 
at the origin leads to spurious transmission 

2^- { kh-ym (61) 

\J 1 — Jkh~yJ\2 + yj\ - (fc*+) 2 /12 

which is valid in the range of resolution in which the weighted-average 
formulation represents propagation (along grid lines). 

HO Transmission for the higher-order representation with the parametei 
defined in (15) is 


„ ^1 -(fc*-)V6 

, _ -(1-(H-)V12) (62) 

^/l — (fc/i“) 2 /6 s]\ - (M+)7 6 

(l-(« 1 -)V12) + (l-(«*+)V12) 

which is valid in the range of resolution in which the higher-order for- 
mulation represents propagation (along grid lines). 

EP The case that is non-dispersive in one dimension (13) may be treated 
similarly. In this case the transmission is 


4>o — 


sin (kh~)/h~ 

2 + cos (kh ~ ) 

sin (kh~)/h~ sin (k,h+)/h+ 

2 + cos(fc/i“) 2 + cos(fc/*+) 


(63) 


which is valid in the range of resolution in which the higher-order for- 
mulation represents propagation (along grid lines). 


Spurious transmission of the various formulations at different wave resolu- 
tions is plotted in Fig. 8. In general the sensitivity to transition in grid size is 
higher for coarser grids. The weighted-average representation is significantly 
superior to the pointwise scheme on non-uniform grids. The higher-order 
and exact-phase formulations offer further improvement. 





5.2 Interface of physical properties 


Discontinuities in physical properties give rise to wave reflection and trans- 
mission. The relative amplitudes of the reflected and transmitted waves 
depend on the ratio of wave numbers, which defines the character of the 
discontinuity. The numerical representation of these phenomena by finite 
element methods was studied in [18]. 

Consider a generalization of the previous configuration in which a discon- 
tinuity in material properties as well as a jump in grid size may occur at the 
origin, so that k~ is the wave number to the left of the origin and k + — to 
its right. An incident plane wave of unit magnitude traveling in the positive 
direction exp(zA:’“x) for x < 0 generates reflected and transmitted waves, so 


that 


4 > = 


exp (ik x) + (<f>(0) — 1)/ exp(ik+x), x<0 
<j>(0) exp(z/c + x), x > 0 


(64) 


where 


m = 


2 k' 


(65) 


k~ + k + 

The discrete solution is again (56) where the numerical wave numbers are 
k h± = k h (k ± h ± ). 


PT The pointwise representation (5) of the solution at the origin yields 

n 2 (fa - (1 - (k + h + ) 2 /2)<t> 0 <f >- 1 - (1 - ( Arfr-) 2 /2)<fto \ 

° " h+ + h~ ( h+ + h- ) 

* ( ( “"(****+> I ^sm(knn\ 

- h+ + h-\\ h+ h- ) 90 h- ) 

where, again, the dispersion relation for the pointwise representation 
(35) is employed. Thus 

2*- v /TuFFjv4 

<j > o — / / (^ ' ) 

k-yfl - (JPh-yi 4 + k+yj 1 - (k+h +) 2 / 4 
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WA Similarly, the weighted-average representation (9) of the solution (56) 
at the origin leads to 


_ 2k-yjl - (k-h-y /12 

k- yj\ - {k-h-) 2 l 12 + k+ y/l - (jfc+A+) 2 /12 


( 68 ) 


HO Transmission for the higher-order representation with the parameter 
defined in (15) is 


2t 


$ o — 


\A - (fc-A-)V6 
(1 -(k-h-f! 12 ) 


yi - U-t-) 2 / 6 yi - (*+A+)V6 

(1 - (k-h-) 2 / 12) + (1 - (A;+/i+) 2 /12) 


(69) 


EP The case that is non-dispersive in one dimension (13) may be treated 
similarly. In this case the transmission is 




— 


_sin(£/i )/(& A ) 
2 + cos(fch”) 


sin(fc/i )/ (A; h ) ^ + sin(fc/* + )/(A; + /i + ) 

2 + cos(kh+) 


(70) 


2 + cos (kh~) 


Physical transmission depends on the ratio of the wave numbers. Nu- 
merical solutions depend on this parameter and on the ratio of resolutions. 
To find out which of the two parameters significantly effects the numerical 
error in transmission consider the transmission error as a function of ratio of 
resolutions for 6 grid points per wavelength to the left of the origin. This is 
plotted for a ratio of the wave numbers equal to unity in Fig. 8 (top). Increas- 
ing the ratio by one order of magnitude and by two yields the behavior shown 
in Fig. 9. The difference between these plots is not significant indicating that 
the error depends primarily on the ratio of resolutions. All the representa- 
tions have the property that <£ 0 = <£(0) if k+h+ = k~h~ . Again, superior 
performance of the weighted-average representation and its enhancements is 
evident. 
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6 Local Truncation Error Analysis 

The local truncation error is the residual left by substituting the exact solu- 
tion in the discrete representation. In the following, sufficient differentiability 
is assumed for all functions involved. 


6.1 Uniform grids 

Consider the one-dimensional constant-coefficient case on a uniform grid. For 
the pointwise representation (3) 

r = ~ + + e<t ,( z ,) + /<*,) 

= n *,) + + <roo) + + n*>) 

= ^r(ii) + 0(ft 4 ) (71) 

where primes and superior Roman numerals indicate differentiation by the 
argument. The second line is obtained by Taylor’s formula, where > 
x~ > Xj and Xj > > Xj+ 1 , and the third line, which follows from the fact 

that <f> satisfies the Helmholtz equation, indicates consistency. The pointwise 
scheme is thus second-order accurate. 

The weighted-average representation (6) is similar on uniform grids 

_ _ 1 ) - 2 <f>(xj) + <j>(xj. i) _ j2 (f>{x j+ i) + 4<ft(xj) + , 

+k + 

f( x j+i) + ^/( x j) + /( x j-i) 

= ~ + + + 

6 

= + 0(h<) (72) 

The weighted-average scheme is also second-order accurate. This order of 
accuracy is retained in the case of variable coefficients (7). 
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For the improved representations (14) 


<£(zj +1 ) - 2<j>(xj) + ,2 <ftpj+ 1) + 4^(gj) + 4>{?j- 1) 


/3h 2 


+ k 


+ 


f( x j+ 1 / 2 ) + /pi) + /pj- 1 / 2 ) 


3 


= 1 + 


/ 3(kh) 2 \ <j>(xj+i) - 2<f>(xj) + <f>{xj-i) 


f{ x j) + 


6 ) (3h 2 

f( x j+ 1 / 2 ) — 2 f(xj) + f(xj_ 1 / 2 ) 


+ k 2 (f>(xj ) + 


(73) 


Employing the definition of /3 that leads to the high-order representation (15) 
yields 

T = - r(*l)l 4) + 0(ft s ) (74) 

justifying its name as a higher-order scheme. Note that if the source terms 
were not represented appropriately there would be second-order terms in the 
truncation error. The scheme that is dispersion-free in one dimension (13) 
has a truncation error 


t = ^ (fc2w “ (xj) + < 75 ) 

If the fourth derivative of the source vanishes the method becomes six-order 
accurate. Furthermore, the truncation error is zero when all the derivatives 
of the source from fourth order and higher vanish. 


6.2 Non-uniform grids 

In analyzing method performance on non-uniform grids a change of variables 
from physical space to computational space is often considered, so that the 
grid is uniform in the latter [19]. The order of accuracy of some methods 
on non-uniform grids may drop in physical space. Nevertheless, in computa- 
tional space it remains unchanged from the order in the uniform case. To a 
certain extent grid stretching should reflect variation of the solution. In this 
case accuracy in computational space is representative of the situation. In 
practice, however, grid variation is determined by geometric considerations 
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as well. The following results are thus presented in physical space, which 
describes the more general case. 

For the pointwise representation (4) 


r = 


f <j)(x j+ i) - <f>(xj) <t>{xj) — /h + + h 2 


V 

<t>"{Xj) + 


h+ 


h+ - h- 


■r(xi) + 


(k + ) 2 -h + k- + (h-y 
12 


+ k <f>(xj) + f(xj) 


r ixj) + 


60 (h+ + h~) 


(( h + ) 4 r(x + ) - (h-) 4 p(x-)) + k 2 <t>( Xj ) + f( Xj ) 


h + - h- 

3 


<j>'"( x j) + 


(h + ) 2 - h+h~ +(h-) 2 

12 


■r(xj) + o(h 3 ) 


(76) 


The pointwise scheme is indeed second-order accurate on uniform grids, but 
may drop to first order in the non-uniform case. 

Whether the scheme actually drops to first order or not depends on the 
degree of grid stretching. If 


h + - h~ = 0(h p+ '), p > 0 (77) 


the stretching is called algebraic [19]. With algebraic stretching the pointwise 
scheme retains second-order accuracy. Otherwise the accuracy drops to first 
order. 

In contrast, for the weighted-average representation (8) 


( 4>{x } +\) - 4>( X j) - <f>(xj-i)\ / h+ + h~ 


V ^ 

k 2 ( h+ 


3 \h+ + h 
1 / h + 


h~ 


— (^(xj+i) + 2<jl>(£j)) + 


h+ + h- 


+ 


(2 <j>( Xj ) + ^(xj_,)) + 


3 \ h+ + h 


r(/( x j+i) + 2 f( x j)) + 


h+ + h 


z( 2 f(Xj) +f(Xj-i)) 


2 


h+ + h- 


f(xj) + 


1 + 


(kh + ) 2 \ <t>(x }+ 1 ) - <j)(xj) 


6 


h+ 


1 + 


1 


(M ) 2 \ ^(gj-i) - 4>(xj y 
6 ) h- 


~b k 2 <j)(xj) + 


3 (h+ + h~) 


rr (A + (/(x J+1 ) - /(i,)) + h-(f( Xj ) - /(*,.,))) 
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(78) 


= ^ + (/, - ^o‘"(.c/) + 0(h 3 ) 

first-order terms cancel out even on non-uniform grids so that second-order 
accuracy is retained in any case. Again, these results apply to the case of 
variable coefficients as well. 

For the improved representations 


j{ fi*k^+W wxi *' )+2Hzi)) 

+ & +h++l3-h- mXl) + + 


P + h + 


3 \P + h+ + f3~h~ 
0-h 


— (2f(x J + i/ 2 ) + /(^j)) 


+ ^TrF </(l,) + ' 2/(Xj - ,/2)) 


p+h+ + p-h- 


k 2 fi + h + h + \ <f> (x J+1 ) - <f>(xj) 


k 2 fj h h \ ^(gj-i) - 2 


+ k 2 4>(xj ) + 


/( x j) + 


3(P + h + + p~h~) 


[p + h + {f{x J+l/2 ) - f{ Xj )) 


+ P h (f(Xj) - f(Xj-i/2 


For the definition of p that leads to the high-order representation (15) 

t- = ~ ' r)(( 9o )2 + ( — (*V"(*j) - f"MI‘ i) + 

(t+)* ~(A+) 3 A + (fe + ) 2 (A f-h+jh ) 3 + (A ) _ /••( Xj)/4) + 

240 

0{h 5 ) (80) 

the truncation error is third-order accurate on non-uniform grids (and, of 
course, fourth order in the uniform case). 
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7 Conclusions 


In this work finite difference methods for solving problems of time-harmonic 
acoustics are developed and analyzed. The well-known pointwise represen- 
tation, eqs. (3) and (19) in one and two dimensions, respectively, is second- 
order accurate on uniform grids. However, accuracy may drop to first order 
in the non-uniform case (4) unless sufficiently smooth grid stretching (77) 
is employed. In multi-dimensional configurations the representation actually 
improves the less aligned the propagation directions are *with respect to the 
grid. 

A weighted-average representation, eqs. (6) and (21) in one and two di- 
mensions, respectively, has the same asymptotic behavior on uniform grids, 
but is less sensitive to low wave resolution and, more importantly, to di- 
rection of propagation and transition in wave resolution (including material 
interfaces). Performance in multi-dimensional configurations again improves 
for propagation directions that are not aligned with the grid. In general, 
anisotropy in numerical representation is reduced with increased wave reso- 
lution. At lower resolution the weighted-average representation (21) is much 
less anisotropic than the standard pointwise representation (19). Second- 
order accuracy is retained on any non-uniform grid (8) at virtually no in- 
crease in computational cost. These results hold for variable coefficients as 
well. 

Superior performance is attained by basing the schemes on a generalized 
definition of the derivative (10) which incorporates a resolution-dependent 
parameter. Improved schemes with higher-order accuracy are designed by 
appropriate definition of the parameter (15), reducing spurious dispersion 
and reflection. Defining the parameter for schemes which are, in some cases, 
dispersion-free (13) leads to the same asymptotic behavior with improved 
coarse grid accuracy. Source terms must be represented accordingly (14) so 
as not to degrade the higher-order accuracy. These methods are, in general, 
fourth-order accurate on uniform grids and third order in the non-uniform 
case. The performance of these schemes in multi-dimensional configurations 
is superior for any direction of propagation. Their performance improves 
as propagation directions become aligned with the grid. In principle, grids 
should thus be aligned with directions of propagation to the extent possible, 
further enhancing the performance of these methods. 

Schemes that exhibit higher-order behavior on uniform grids in all di- 
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rections in two-dimensional configurations are derived on the basis of Pade 
approximation and its generalization (22). The dispersion of these methods 
(as well as their spurious reflection and transmission) along grid lines is iden- 
tical to that of the higher-order method based on the generalized definition 
of the derivative. Dispersion along grid diagonals is minimized by employing 
7 = 2/5 which leads to (27). These methods are by far less anisotropic than 
all other schemes. The value of 7 = 14/5 leads to the stencil with the lowest 
degree of anisotropy (28). 

In general, wave resolution (kh) should be kept as even as possible through- 
out the grid to minimize spurious reflection and transmission. Sensitivity to 
these phenomena is greater on relatively coarse grids. 
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